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Abstract. We present a method that assesses the theoretical detection limit of a 
Bayesian Markov chain Monte Carlo search for a periodic gravitational wave signal 
emitted by a neutron star. Inverse probability yields an upper limit estimate for the 
strength when a signal could not be detected in an observed data set. The proposed 
method is based on Bayesian model comparison that automatically quantifies Occam's 
Razor. It limits the complexity of a model by favoring the most parsimonious model 
that explains the data. By comparing the model with a signal from a pulsar to the null 
model that assumes solely noise, we derive the detection probability and an estimate 
for the upper limit that a search, for example, for a narrow-band emission for SN1987a, 
might yield on data at the sensitivity of LIGO data for an observation time of one year. 

PACS numbers: 04.80.Nn, 02.70.Uu. 
1. Introduction 

Several mechanisms have been proposed that would cause rapidly rotating neutron stars 
to emit quasi-periodic gravitational waves [H [2j. Interferometric gravitational wave 
detectors that are now operating in numerous locations around the world [H HI [5j [6] now 
allow for their verification and much work has gone into the development of dedicated 
search algorithms for these signals. Radio observations can provide the sky location, 
rotation frequency and spin-down rate of known pulsars. The frequency of the reported 
remnant of SN1987a for example is not known accurately [7j but Markov Chain Monte 
Carlo (MCMC) methods [HI El El EH EE] are able to search a range of frequencies (and 
other physical parameters) in a reasonable time. 

As in previous studies [T3J [H] the signal under consideration is one that is expected 
from a non-precessing triaxial neutron star. The gravitational wave signal from such an 
object is at / = 2/ r twice its rotation frequency f r , and we characterize the amplitudes 
of each polarization with overall strain factor, ho. The measured gravitational wave 
signal will also depend on the antenna patterns of the detector for the 'cross' and 'plus' 
polarizations, F x>+ (t; ip, a, 5), giving a signal s(t) =F + (t; . . .)/i (l+cos 2 t)cos $(£; . . .)/2+ 
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F x (t; . . .)ho cos i sin $(£; . . .) where i is the inclination angle. The antenna pattern of the 
detector depends on time t, the polarization angle ip and its location determined by right 
ascension a and declination angle 5. The location is assumed to be known from, for 
example, radio observations. A simple slowdown model [15] provides the phase evolution 
of the signal as 

$(t; n, f„ f s ) = O + 2tt [f a (T (ajS) - T ) + / S (T (Q)<5) - T ) 2 /2] , (1) 

where 

T {a , S ) =t + St = t + ^ + AT (2) 

is the time of arrival of the signal at the solar system barycenter when t is the time at 
the detector. Here, 4>o is the phase of the signal at a fiducial time T , r is the position 
of the detector with respect to the solar system barycenter, n is a unit vector in the 
direction of the neutron star (depending on a and 5), c is the speed of light and AT 
contains the relativistic corrections to the arrival time [16] . 

If f s , f s , and n are known from radio observations, for instance, the signal can 
be heterodyned by multiplying the data by exp[— n, f s , f s )], low-pass filtered and 
resampled, so that the only time varying quantity remaining is the antenna pattern 
of the interferometer. The reference sky location is also needed for the heterodyning 
process prior to the MCMC simulation. We are left with a simple model with four 
unknown parameters ho, ip, <fio, and l. If there is an uncertainty in the frequency and 
frequency derivative two additional parameters come into play, the differences between 
the signal and heterodyne frequency and frequency derivatives, Af and Af. The unit 
vector n points to the right ascension a and declination S of the purported neutron star. 

A detailed description of the heterodyning procedure is presented elsewhere [TBI HI] . 
The model of the heterodyned signal of a pulsar has form [H] 

y(t k ; a) = F + (t k ; ^ a, 5)h (l + cos 2 t yA*fe<*,W,A/)/ 4 

- tF x (t k ; V>, a, 5)h cos W,Afy 2j (3) 

where t k is the time of the k th bin and a = (ho, cost, (fio, if), Af, Af) is a vector of 
the unknown parameters. A$(i; a, 5, Af, Af) represents the residual phase evolution of 
the signal, equaling O + 2tt[A/(T (Qi(5) - T ) + A/(T (a>a) - T ) 2 /2], where T M (Eq. P) 
depends on the known sky location of the pulsar. Note, that the gravitational wave 
oscillates at twice the rotation frequency of the pulsar's rotation frequency. Therefore, 
the frequency in Eq. [3] refers to the gravitational wave frequency. The objective is to fit 
this model to the data B k = y{t k ; a) +e k , where e k is assumed to be normally distributed 
noise with a mean of zero and known variance a\. Assuming statistical independence 
of the binned data points, B k , the joint likelihood that these data d = {B k } arise from 
a model with a certain parameter vector a is [H] 

p{d\a) ex l[exp \- \(B k - y{t k ; a))/a k \ 2 } /2 = exp [-x 2 (a)/2j , (4) 



where 



X 2 (a) = Y.\B k -y(t k ;a)\ 2 /ai (5) 

k 
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In order to draw any inference on the unknown parameter vector a the posterior 
probability of a given d is needed, which can be obtained from the likelihood via an 
application of Bayes' theorem. The unnormalized posterior density p(a\d) oc p(a)p(d\a) 
is the product of the prior density of a, p(a), and the joint likelihood, p(d\a). In this 
study uniform priors distributions are used with prior ranges [0, 2ir], [— ir/4, 7r/4] and 
[— 1, 1] for the angle parameters (fi , ip and cost respectively. 

For ho, a uniform prior is specified with boundary [0, 10~ 20 ]. For the frequency and 
spin down uncertainty, suitable uniform priors are used with ranges of [— y^j, y^j] Hz and 
[— 10~ 9 , 10~ 9 ] Hzs -1 for Af and Af, respectively, as applied in [10]. The normalized 
posterior density p(a\d) = p(a)p(d\a) / p(d) cannot be evaluated analytically, therefore 
Monte Carlo methods are used here to explore p(a\d), as described in [TO] . 

When the signal-to-noise ratio (SNR) and hence the signal's evidence declines, 
it becomes increasingly difficult to sample efficiently from the posterior distribution 
using MCMC The major problem lies in the frequency parameters Af and Af. Long 
integration periods yield narrow posterior modes and when the SNR is small, their 
occurrence is also negligible with most of the posterior probability mass spread over 
the entire parameter space determined by the prior distribution. The sampling process 
of an MCMC sampler becomes inefficient in covering that part of the parameter space 
where the signal is concentrated. The question that will be addressed in this paper 
is the threshold of the SNR for which MCMC sampling becomes ineffective and below 
which no signal parameters can be retrieved. 

2. The detection of weak signals 

The presence of a signal within the data can be assessed by a formal Bayesian model 
comparison of the model the contains a signal with the null model that contains no signal. 
Bayes factors could be applied but they require a properly converged MCMC output. 
Without the need of MCMC samples, this paper aims to give theoretical detection 
probabilities dependent on signal-to-noise ratios. 

2.1. Derivation of a theoretical detection probability 

For the Bayesian Information Criterion (BIC), also called the Schwarz criterion, there 
is no particular need for the MCMC output samples. The BIC is defined [17] as 

BIC = —2 log(maximum likelihood) + V, (6) 

where the penalty term V = dlogn brings in the number of d = 6 independent 
parameters that describe the model, and the number n of data samples. The penalty 
term penalizes the number of parameters in a model in order give preference to simpler 
models and meet the principle of Occam's Razor. 

The objective is to derive a theoretical limit for the detection of a signal within a 
data set observed during a determined observation period at a certain noise level. This 
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section is dedicated to find a distribution of the BIC depending on the noise, conditioned 
on the parameters of a potential pulsar. 

The observation period is a vector OP = (ti,...,t n )' of n time points tk with 
k G {1, . . . ,n} during which the data has been collected starting from t star t and ending 
at tend- The noise vector is a vector cr = (ci, . . . , a n )' for the n data bins. Given the true 
parameter vector of the pulsar from which the signal arises, the full information needed 
for a detection is determined by the vector a* = (h^, cos t*, ip*, a*, 5*, A/*, A/*, cr, OP)'. 
Although some parameters like the sky location are expected to be known, they are 
essential factors for the detection probability in connection with the observation period 
and the noise. These are essential parts of the parameter vector as the detection depends 
significantly on them. 

A signal detection depends on the actual evidence of the model that assumes the 
presence of a signal from a pulsar within an arbitrary data set when compared to the 
null model of mere noise. Each potential data set under consideration is based on the 
true parameters of a potential pulsar. Therefore each model comparison is conditioned 
on a data set that is conditioned on the parameter vector a*. This fact can be used 
to obtain, for large sample sizes, an approximation for the maximum likelihood value 
since the maximum likelihood estimate (MLE) is asymptotically consistent and efficient 
under certain regularity conditions that are generally satisfied [15] . Thus the estimates 
converge to the true values for large samples sizes. The sample sizes that we expect are 
in fact in the range of tens of thousands. 

A potential data set d* from a pulsar, based on a true parameter vector a* is 
modeled by : = y(tk, a*) + e& with noise vector e^. Due to the fact that is 
conditioned on a*, an approximate maximum log-likelihood under model Aii is 



This term comprises the sum of the squared residuals as the model is fitted by the 
true parameter vector. On the other hand, under model Aio that encompasses no 
parameters, the log-likelihood has a constant value and therefore its maximum is 



where the summation term contains the true and given parameter vector of the signal. 
It is clear that logMLd», a »,Xi > logMLd^a.^o^' 2 *- As a resu h of this, naturally model 
Aii has to be preferred at all times. This, however, does not take into account the 
penalty term that comes into play due to the principle of Occam's razor. Equality 
of Eq. [8] and [7| can only be achieved for a zero amplitude Hq in parameter vector a*. 
But how large do we have to choose this amplitude, also considering other influential 
parameters, in order to justify model Ai\ with its many more parameters? This is the 
essential idea behind this model comparison approach and the penalty terms play a key 
role in it. 
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We aim to compare model M. and .Mi conditioned on the data set d*, conditioned 
on a potential pulsar characterized by the true parameter vector a*. By substituting 
Eq. (ED and Eq. (ED into Eq. ©, we obtain 

BICd, )0 , iM) = -2 log ML d ^ atMo (9) 
as model .Mo has d = parameters and 

BKU, „ Ml = -2 log MLd,^^ + P. (10) 
With respect to .Mo and .Mi, a probability for model .Mi can be derived by 

p{Mi\d.,aJ = (l + e ABIC *.a,/2-io gP (^i)+io gP (^o))- 1 (n) 

Here, p(-M ) an d p(A^i) are prior probabilities for A4 and .Mi respectively The 
interested reader is referred [19] for a more detailed derivation. We will address different 
prior scenarios later but for now, we choose equal probabilities p(M. Q ) = p(M.i) = 0.5 
for the models as a natural choice when there is no prior information about the 
possible existence of a signal. This yields p(M.i\d*, a*) = (l + e ABICd *. a */ 2 ^ where 
ABICd^a, := BlCd t ,a t ,Mi ~ BICd»,a»,A4 - ^ represents the probability that the data 
from a potential pulsar with given parameter vector a* is better modeled by Aii (a 
signal) rather than .Mo (no signal). In other words it is the probability for the existence 
of a signal in the data that is emitted by a pulsar with parameter vector a*. It is 
merely the difference of the two BIC values under consideration that is responsible for a 
signal detection. A difference of zero for example would yield a 50% probability for both 
models. A probability conditioned on data d* from the vector a t , can be expressed as 

p(Mi\a*) = E \p(Mi\d*, a*)|a*] = E 

There is no simple way to solve this expression analytically and although feasible, 
a Monte Carlo sampling process would be lengthly. From a physical perspective, phase 
0o and the frequency parameters Af, Af should have no impact on the actual signal 
detection as the SNR mainly depends on the amplitude Hq, inclination cos i*, noise <x, 
and observation time OP. To a smaller extent the SNR is also influenced by the course 
of the antenna pattern over the observation time OP with parameters ip*, a*, and 5*. 
We assume the sky location to be known and condition on a* and 5*. 

The probability p(M.\ \ a*) is determined by the distribution of ABIC^a,. Thus the 
characteristics of ABICd,^, will be derived below. By using equations Eq. ©, Eq. ©, 
Eq. ©, Eq. Jffi) we obtain 

ABic*,^ ~EM 2 K 2 + ^-E + e *l 2 f°k (is) 

k k 

In [14|, white Gaussian noise e^re, efc jim ~ iV(0,crj?) is assumed where the a\ 
are estimated for each bin from the noise floor in a 4 Hz band of data around the 
signal frequency. By substituting y(t k ; a*) of Eq. [3] and defining some abbreviations, 
F +iX (t k ;^*,a*,5*) := F+' x , e ^Ht^*,5*ArAn := e ^ k = cos (A$ fc ) + j s in(A$ fc ), 



l + e 



ABIC 



/2\-l 



a, 



(12) 
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\h^(l + cos 2 l*) =: A + , and ^h^cost* =: A x we can rewrite Eq. (ITS]) as 



-2£(K 

k 



2 



E[(^ + )7^ + (^) 2 M 5 



^ + cos(A$ fc ) + A x F k x sin(A$ fc ) /<j fc e Km /a k 



^+sin(A$ fe ) 



A X F, 



k cos(A$ fc ) /a k ) e kM /a k . 



(14) 



The quadratic noise terms cancel out and we are left with normally distributed terms. 
Given a pulsar with parameter vector a*, the ABIC,^^ is thus normally distributed. 
The terms that contain the phase evolution canceled out as well and Eq. ( fl~4"l) is thus 
independent of the parameters (fi , Af, and Af. With e^re, e^im ~ N(0, o~ k ) we have 
E^re/c/c) = E(efc 5 i m /o"fc) = and the expected value of Eq. ( fl~4"l) has the form 



/V := E(ABIC d „ a J = V - \(A + F+r + (A X F, 



tx\2 



(15) 



Eq. (Tl5|) yet allows some insight as it tells us that for a given arbitrary parameter 
vector a*, model A4q would be preferred over Mi, if fi at > 0. Given a parameter vector 
a*, the variance of Eq. (fl4l) is 

< = Var(ABIC d „ a J = A(V - /i a J. (16) 

Both expressions Eq. ([15!) and Eq. (fl6|) only depend on the five parameters h$, 
cos 6*, -0*, a*, and 5*. The parameters ip*, a*, and 5* only enter in the plus and 
cross polarization terms F k and F k x of the antenna pattern which depends on the 
orientation sweep of the interferometer towards the pulsar and the polarization angle of 
the gravitational wave that it emits. 

We are left with the random variable ABIC|a* ~ N(/i aji ,cr 2 J that depends on 
five parameters of the pulsar plus noise cr and observation period OP. If we assume 
constant noise a over time, we can combine h* and a to a more handy SNR h^/a 
parameter. We define a new vector a. = {hi/ a, cos l* , i/;*, a*, 5*, OP)' with observation 
period OP = (ti, . . . ,£„)'. Explicitly, the difference in the BIC values with respect to 
models Mo and Mi, for arbitrary data sets, conditioned on a. follow the distribution 
ABIC|a. ~ N(>bic .a, ■ ^BIC a.) "With 

2 



/^BIC.c 




4 



1 + cos 2 i* 



COS L 



E (17) 



and 



0~ t 



'BIC,a. - 4 (P - /iBIC,a.) • 

Using these information, Monte Carlo methods can be used to estimate Eq. (FEZ]) . 

As an example, we consider a data set that we can expect was taken over one year at 
the three LIGO interferometers Hanford (4km, 2km) and Livingston (4km) with three 
different noise levels at the three interferometers. A sensible heterodyning frequency 
for the SN1987a remnant is / s = 2f r = 935Hz [7]. For the purpose of illustrating an 
example we will assume noise levels that are likely to be close to LIGO's S5 values at 
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the frequency in question. We therefore assume noise levels 8 x 10~ 24 (Hanford 4km), 
1.5 x 10 -23 (Hanford 2km), and 9 x 10~ 24 (Livingston 4km) at that frequency and an 
observation period OP of one year of S5 data that would be heterodyned to a potential 
source at a* = 5 h 35 m 28.03 s and 5* = -69° 16' 11.79" (SN1987a) with 525600 bins at one 
sample per minute. The data are analyzed for each interferometer separately and also 
combined by the sum of the log-likelihoods, as we assume independence. The parameter 
vector encompasses a, = (h*/cr, cos l\ ip*, a* = 5 h 35 m 28.03 s , 5* = -69° 16' 11.79", OP) in 
which the values of h^/a and cos l* and ip* are unknown. In order to derive a probability 
conditioned on /iq/ct, we need to marginalize p(M-i\a 9 ) over cost* and ip and obtain 
p(Mi\h*/a,a*,5*,OP) = J p(Mi\a,)dp COSL (cos L*)dp^(ip*). 

Fig. [1] displays the probability of a signal detection as a function of the amplitude. 
Two different prior probabilities on the signal existence are chosen. The natural choice 
is p(A4i) = 0.5 when there is no information available. However, we know that we 
focus only a l/60Hz band and the probability of an existence needs to be split on the 
frequency bands in which we expect a signal. In addition, we do not know whether there 
is a neutron star at all which lowers the probability further. For this reason, we chose a 
rather arbitrary and extremely small probability of p(Aii) = 10 -9 in order to asses the 
impact of that prior probability. We obtain the graph shown in Fig. [IJ 



Figure 1. Signal detection probability 
for the three different interferometers 
and two different prior probabilities 
for model A^i as a function of the 
amplitude Hq for one year of S5 data. 
The curves of two prior probabilities 
p(Mi) = 0.5 (dashed lines) and 
p{M.i) — 10 -9 (solid lines) are shown. 



A larger amplitude h$ is required for a successful detection when we doubt the 
existence of a signal. Hence, the data must speak more clearly for a signal in order to 
overcome the low prior probability but since the observation period of one year is rather 
long, the effect of the prior probability is fairly small. 

All graphs compiled so far are showing a signal detection probability given a 
particular scenario but the question we aim to answer in the next section is how strong 
a signal still can be even if a signal can't be seen. 

2.2. Performance of the Bayesian MCMC search in setting an upper limit 
using S5 data 

The upper limit estimate for a Bayesian MCMC search involves testing the hypothesis 
/iq < UL vs. /iq > UL under the assumption Ai that there is no signal in 
the data. The derivation of the probability p(h$ < UL|A^o) will shed light on 
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this matter. We condition on noise, observation period, and location and after 
integrating over the prior distributions of cos t* and ip*, p(A4o\hQ,OP,er,a*,5*) : = 
/ / p(-M o I a , )dp COS L (cos i*) dp^p (if)*), we obtain 

v(h* < ULIM OP a a* 5*) - C p(Mq\K, OP, cr, a*, 6*) P (h*)dh* 
p(h < VL\Mo, OP, a, a , 5) - ^ p{M ^ 0Pj - a%5 . )m)dK ^ 

In order to derive Eq. ( fT9l) we need to find a suitable prior for p{h* Q ). One choice could 
be to put a uniform prior on Hq with large boundary [0, 10~ 20 ]. The upper boundary of 
the prior range has negligible impact on the results of Eq. ( 1T91) as long as this boundary 
is significantly larger then the upper limit estimate. Fig. [2] displays Eq. (|T9|) for two 
different prior probabilities on whether we expect a signal at SN1987a. 




Figure 2. Estimated sensitivity of the Bayesian method described in this paper, 
assuming 1 year of data with the typical noise level of the LIGO interferometers during 
their S5 run. The model prior probabilities are p(M.i) — 0.5 (left) and p(Aii) = 1CP 9 
(right). The prior for ho is ho ~ Unif(0, 10 -20 ). The assumed noise levels are 
cr Hl = 8 x 1CT 24 , cth 2 = 1-5 x 1(T 23 , and cr Ll = 9 x 10" 24 . 



Since we focus our search on a possible pulsar in SN1987a, we can tailor a prior 
distribution for ho as we know the age of SN1987a and its distance. In [20] it is assumed 
that a newly formed neutron star spins at high rate and gravitational radiation slows it 
down. Two different prior scenarios are conceived here. According to [20], it is 

h(f) = r- 1 V / (5G^)/(8c 3 r gw (/)), (20) 

where T gw (f) is the time for the gravitational wave frequency to drift down to frequency 
/ from its original spin rate. In case of SN1987a it is 20 years. Here, G is Newton's 
constant, c the speed of light, r the distance to the neutron star, and I zz the principal 
moment of inertia about the rotation axis. In order to derive a prior distribution for ho 
we need to determine prior distributions for r and I zz . We assume the distance estimated 
in [21] with 50.9 ± 1.8kpc with r ~ N(50.9, 1.8 2 )kpc for accounting the uncertainty in 
the distance. For the moment of inertia, we choose a uniform prior within the range 
[10 38 , 3 x 10 38 ] as applied in [22]. These considerations yield a prior for ho as shown 
later in Fig. 

A totally different approach for obtaining a prior distribution for h is by [13] 

ho = An 2 GI zz fe/(c 4 r) (21) 
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for a general pulsar expected at SN1987a. Here, / is the pulsar's rotation frequency, 
and e its ellipticity. We heterodyne to a frequency of / s = 2/ r = 935Hz, and assume 
the gravitational wave frequency to have this value within a 1/ 60Hz frequency band for 
a particular search. An uncertainty beyond this needs to be accounted for in the prior 
p(A4i) for the existence of the signal within the l/60Hz band around / s because the 
signal is not seen outside that band after the heterodyning process. 

We use the same uniform prior for I zz as above but we have to find a suitable prior 
for the ellipticity. In [23], the ellipticity is assumed to have an exponential distribution 
(maximum entropy prior) with cut-off at a maximum ellipticity threshold. Although in 
[23] more pessimistic mean and maximum values are used, our choice is more optimistic 
in order to account for the fact that we know that a possible neutron star in SN1987a 
is very young. We choose a cut-off according to [20J at e max ~ 9 x 1CT 5 based on the 
idea of a hybrid neutron star with a mixed quark and baryon core and a normal neutron 
star in the outer part. For the mean of the exponential prior distribution we use an 
optimistic choice of e mean = 5 x 1CT 5 . Both prior distributions for ho as discussed above 
are displayed in Fig. [3] along with their resulting upper limit estimates. 
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Figure 3. Two different prior distributions (ieft column) for ho for a possible neutron 
star in SN1987a and the corresponding curves for the Bayesian MCMC upper limit 
estimates (right column). The upper row corresponds to a prior subject to Eq. (I20p 
whereas the lower row is based on Eq. (j!?Tj) . For the model selection, the prior 
probability is chosen to be p(Mi) = 1CT 9 . 

The use of such priors changes the results for upper limit estimates compared to 
those in Fig. [2] (which were based on a uniform prior). In essence, a uniform prior 
on the amplitude recovers the detection ability of an interferometer. For example, in 
case of combined data sets, an upper limit estimate based on a uniform prior requires 
an amplitude of at least 6.2 x 10~ 25 . The use of prior distributions based on Eq. ( |20l ) 
and Eq. (!2T|) . however, only have 0.001% and 11.4% probability mass above that limit, 
respectively. This inevitably yields values for the upper limits estimates dominated by 
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the prior of ho. This is obvious especially in case of the prior based on Eq. (1201) and can 
be seen in Fig. [3j 

3. Conclusions 

The Bayesian MCMC methods work well when the SNR is sufficiently large but they 
struggle when the signal is too weak and the parameters that affect the phase evolution 
are not known. The fact that we integrate over very long observation periods requires 
an almost exact match of the phase evolution and almost all mass of the posterior 
distribution is highly concentrated around one point in the parameter space when the 
SNR is large. Finding this posterior peak with Bayesian MCMC methods is time 
consuming but once found, the sampling process is easy and efficient. With decreasing 
SNR, however, the sampler is forced to also sample from other areas of the parameter 
space determined by the prior. This requires multiple retrievals of the narrow peak and 
it requires extremely long runs to gain insight into the actual shape of the posterior 
distribution. The sampling speed depends on observation length and number of Markov 
chains involved when using parallel tempering. For one year of data, each single chain 
samples about 150 000 samples per week and chain on a 2.8 GHz machine. At low SNRs 
at least 10 chains are needed [19]. When no sensible inference can be drawn from an 
MCMC output if no frequency parameters can be retrieved. In those cases, the method 
derived here, based on model comparison, provides an excellent means for estimating an 
upper limit for the amplitude of a signal when using Bayesian MCMC methods, given 
the observation period and noise. In practice this method could be used to estimate 
the sensitivity of the Bayesian MCMC method on actual S5 data and in particular 
the noise. For long observation periods, the impact of prior information about the 
presence of a signal is rather small. The influence of the amplitude's prior only becomes 
significant when the sensitivity, with respect to the obtained data, is too small for the 
expected amplitudes. Consequently, when we expect amplitudes below the detection 
limit then the upper limit estimate is determined mainly by the prior distribution of the 
amplitudes. 
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